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We consider the low-temperature thermal transport properties of the 2D proximity-induced su¬ 
perconducting state formed at the interface between a 3D strong topological insulator (TI) and a 
d-wave superconductor (dSC). This system is a playground for studying massless Dirac fermions, as 
they enter both as quasiparticles of the dSC and as surface states of the TI. For TI surface states 
with a single Dirac point, the four nodes in the interface-state quasiparticle excitation spectrum 
coalesce into a single node as the chemical potential, p, is tuned from above the impurity scattering 
rate (|/r| 3> To) to below (|p| <C Fq). We calculate, via Kubo formula, the universal limit {T —>■ 0) 
thermal conductivity, ko, as a function of p, as it is tuned through this transition. In the large and 
small 1^1 limits, we obtain disorder-independent, closed-form expressions for ko/T. The large-|/r| 
expression is exactly half the value expected for a d-wave superconductor, a demonstration of the 
sense in which the TI surface topological metal is half of an ordinary 2D electron gas. Our numeri¬ 
cal results for intermediate |/r| illustrate the nature of the transition between these limits, which is 
shown to depend on disorder in a well-defined manner. 

PACS numbers: 74.25.fc, 73.20.-r, 74.20.Rp, 74.45.-l-c 


I. INTRODUCTION 

Topological insulatorsip— (TI) represent a novel state 
of quantum matter that comes about due to the com¬ 
bined effects of spin-orbit interactions and time-reversal 
symmetry4p— Though characterized by a bulk band gap, 
they are adiabatically distinct from ordinary insulators 
and support protected gapless surface states. In the case 
of a three-dimensional strong TI, these surface states 
form a novel two-dimensional topological metal with a 
spin-polarized massless Dirac energy spectrum. The the¬ 
oretical prediction and subsequent experimental discov¬ 
ery of TI states in 2D material a^d° (HgTe/CdTe quantum 
wells), 3D material&ldJ^ (BLSbi-a,), and the cleaner, 
simpler, second generation 3D materials^^^— (Bi 2 Se 3 , 
Bi 2 Te 3 , and Sb 2 Te 3 ) has led to great interest in this 
area, exploring both the fundamental physics as well as 
the potential for applications to fault-tolerant topological 
quantum computatio n^^d'^ . 

The proximity of either magnetic materials or super¬ 
conductors to the TI surface can induce an energy gap 
in the topological metal, resulting in even more ex¬ 
otic interface statesii Early on, Fu and Kanei^ con¬ 
sidered the proximity effect at the interface between 
a TI and a conventional s-wave superconductor, ana¬ 
lyzing the proximity-induced superconducting interface 
state and finding that it should support Majorana bound 
statesi^^— at vortices. Subsequent work has expanded 
this analysis in many directions, and has included the 
case of TIs proximity-coupled to unconventional super¬ 
conductors of different pairing symmetries^lvS^. Such 
TI interface state superconductivity has been demon¬ 
strated, experimentally, both for the s-wave cas o^^'^^ and 
for the case of TIs coupled to high-Tc cuprate d-wave 
superconductors^. 

It is this last case, that of the proximity-induced super¬ 


conducting state at the interface of a 3D strong topolog¬ 
ical insulator (TI) and a d-wave superconductor (dSC), 
that is our focus here. For simplicity, we will consider a 
TI with a surface state characterized by a single Dirac 
point at the origin of fc-space, as is seen in the Bi 2 Se 3 
family of material a^^d"^ . We are particularly interested 
in the low energy quasiparticle excitations of this inter¬ 
face state, a system in which massless Dirac fermions 
enter in two different ways, as both the surface states of 
the TI and the quasiparticles of the dSC. For the for¬ 
mer, the TI surface states, the massless Dirac fermions 
are isotropic, a consequence of band structure, and not 
pinned to the Fermi surface, such that one can tune 
through the Dirac point by varying the chemical poten¬ 
tial. They are described by a Dirac equation where the 
gamma matrices live in 2 x 2 spin space. For the latter, 
the dSC quasiparticle states, the massless Dirac fermions 
are anisotropic, their energy spectrum squeezed in k- 
space, and are pinned to the Fermi surface at four nodal 
points. They are described by a Dirac equation where 
the gamma matrices live in 2 x 2 particle-hole (Nambu) 
space. The Tl-dSC interface state that we consider mixes 
both spin and particle-hole space and will have quasipar¬ 
ticle excitations of its own, with features inherited from 
both of the above. 

A useful probe for studying massless Dirac quasipar¬ 
ticles in d-wave superconductors has been low tempera¬ 
ture thermal transport^^— , measurements of which can 
be extrapolated to the particularly simple and interesting 
regime where temperature, T, is small compared to the 
impurity scattering rate, Tq. This is known as the uni¬ 
versal limit because thermal conductivity due to mass¬ 
less Dirac quasiparticles has been shown to be insensi¬ 
tive to disorder in this very low temperature regime42ri^ 
In this paper, we examine the nature of the low energy 
quasiparticle excitations of the Tl-dSC interface state 
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by calculating the universal-limit thermal conductivity, 
KqIT, as a function of chemical potential, /r. Though the 
Hamiltonian for this interface state couples particle to 
hole and spin-up to spin-down, its quasiparticles carry a 
well-defined heat. Thus, thermal transport tracks quasi¬ 
particle transport, and is therefore well-suited to probing 
the excitations of this system. For |/r| <C Tq, it probes the 
single isotropic Dirac node inherited from the TI surface. 
For |/r| ^ Fq, it probes the four anisotropic Dirac nodes 
resulting from proximity-induced d-wave superconductiv¬ 
ity. We study both of these regimes and the transition 
between them as four nodes coalesce into one. 

We begin in Sec. HJby writing down the 4x4 Hamil¬ 
tonian for the proximity-induced interface state, which 
mixes the spin-space Dirac equation of the TI surface 
with the particle-hole-space Dirac equation of the dSC, 
and then solve for the quasiparticle excitation spectrum. 
In Sec. mi we calculate the matrix spectral function, 
derive the thermal current operator, and then use both 
of these to calculate the universal-limit thermal conduc¬ 
tivity tensor via diagrammatic Kubo formula. Closed 
form analytical expressions for kq/T are obtained in 
both the large- 1/i I and small- |p,| limits, both of which 
are discussed in Sec. El Numerical results charting the 
disorder-dependent transition between these two limits 
are presented in Sec. m Conclusions are discussed in 

Sec. ED 

II. PROXIMITY-INDUCED INTERFACE STATE 
A. Hamiltonian 

We consider the proximity-induced superconducting 
state at the interface of a 3D strong topological insu¬ 
lator (TI) and a d-wave superconductor (dSC). For a TI 
like those in the Bi 2 Se 3 family, characterized by surface 
states with a single Dirac point at the F-point of the 
Brillouin zone, the TI surface state is described by the 
Hamiltoniani^ 

Ho = ^il)l{va ■'k- ( 1 ) 

k 

where ipk = (cfe-|-, are electron annihilation opera¬ 

tors, V is the slope of the Dirac cone, fj, is the chemical 
potential, a = ((Ti,(T 2 ) are spin Pauli matrices, and we 
have adopted units where h = 1. Proximity to a dSC 
induces d-wave superconductivity and results in an in¬ 
terface state Hamiltoniani^i^ii^ that is most compactly 
expressed in the following 4x4 Nambu notation 



(2) 

k 


Elk = (vd -k - fi)T 3 + AfcTi 

(3) 

^1= 

(4) 


where the r are particle-hole Pauli matrices that mix 
the ipk and blocks of and the factor of 1/2 
compensates for particle-hole double counting. Here, 
the proximity-induced superconducting order parameter, 
Afc, is of da; 2 _j ^2 symmetry and is taken to be real. (Note 
that in addition to this spin-singlet d-wave term, the form 
of the TI surface Hamiltonian allows for an additional, 
subdominant, spin-triplet (B 2 u) p-wave term to also be 
induced via proximity to a dSCi^ However, as shown 
by Linder et a spin-triplet p-wave pairing am¬ 

plitude in a TI only renormalizes the chemical poten¬ 
tial and never gaps the surface energy spectrum. Thus, 
while its inclusion here would likely result in a quantita¬ 
tive correction to the effect of the singlet term, it is not 
expected to change the essential physics. Thus, for sim¬ 
plicity, we shall defer consideration of the triplet term to 
future work.) Expanding Eq. ([3]) by evaluating the outer 
products of the Pauli matrices yields the 4x4 Hamilto¬ 
nian 


-E 

vk 

Afe 

0 ■ 


vk'^ 

-E 

0 

Afc 

(5) 

Afc 

0 

E 

—vk~ 

0 

Afc 

—vk~^ 

E - 



where ± iky. 

B. Quasiparticle Excitation Spectrum 

The quasiparticle excitation spectrum of the interface 
state is obtained by solving for the (positive) eigenvalues 
of Hk- As shown in Ref. [H for the s-wave case, the 
resulting spectrum is 

Ek = Y/(±u|k| - fj.)^ + Al (6) 

Though the precise functional form of is material- 
dependent, we can proceed, quite generally, as long as 
Afc satisfies two criteria: (1) It has dj; 2 _y '2 symmetry and 
therefore changes sign along the lines ky = Ek^- (2) It 
vanishes faster than linearly with k as k ^ 0. If these 
criteria are met, the quasiparticle spectrum will have the 
following properties. 

For large |/i|, there will be four nodal points in /c-space, 
located at ±kx = iky = where one of the two 

branches in Eq. m goes to zero and quasiparticles can 
be excited for zero energy cost. In the vicinity of each of 
these nodes, 

Ek ~ \/v‘^kl + v\kl (7) 

where v/^ is the slope of A^ at the node and ki and k 2 
define a local coordinate system, centered at each node, 
with the fci-axis perpendicular to the local Fermi surface 
(pointing away from the origin of fc-space) and the fc 2 -axis 
parallel to the local Fermi surface (pointing in the direc¬ 
tion of increasing A^). For energies small compared to 
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the surfaces of constant energy are ellipses, elongated 
parallel to the local Fermi surface for v > (as is typical 
in cuprate superconductors). The presence of disorder 
smears out the nodes, exciting quasiparticles of energy 
less than or on the order of the impurity scattering rate, 
Fq. For T ^ Fo, quasiparticle transport is dominated by 
these disorder-induced quasiparticles which reside within 
ellipses of semi-major axis Fq/ua and semi-minor axis 
Fq/u about each of the four nodes. 

The nodes are distinct for \fi\ :S> Fq, but as \fi\ 
decreases, the inter-node separation decreases, and for 
|/x| <C Fq, the nodes coalesce at the origin of fc-space. As 
long as Afc vanishes fast enough with decreasing k, as 
per condition (2) above, this transition reveals the un¬ 
derlying massless Dirac spectrum of the TI surface state. 
Thus, for \fj,\ -C Fo, 


Ek ~ ulkl (8) 

and the system thereby trades the four anisotropic nodes 
at nonzero k for a single isotropic node at the origin. 
Note that this single node is, however, doubly degenerate, 
as it derives from both of the branches in Eq. ([S]). For 
\^\ and T small compared to Fq, quasiparticle transport 
is dominated by the disorder-induced quasiparticles that 
reside within the circle of radius Fq/u about this isotropic 
node. 


III. TRANSPORT CALCULATION 

Following the approach employed in Refs. and 113, 
we now proceed to calculate the universal-limit quasipar¬ 
ticle thermal conductivity for this system, as a function of 
chemical potential. Key inputs to this calculation are the 
spectral function and thermal current operator, which we 
will calculate first and then utilize in our calculation of 
the thermal conductivity. 


A. Spectral Function 

To obtain the spectral function, we begin by calculat¬ 
ing the Matsubara Green’s function. Working in our 4- 
component Nambu basis, the 4x4 bare Green’s function 
is obtained by inverting the Hamiltonian 

G°(k,fw) = (9) 

where Hk is the 4x4 Hamiltonian from Eq. ([5]). The 
dressed Green’s function is then found via Dyson’s equa¬ 
tion 


where S is the Matsubara self-energy matrix. The re¬ 
tarded Green’s function is then obtained by continuing 
iu! ^ uj + iS 

G^(k,w) = [wl - S^(a;) - Hk]~^ (12) 

where = T,{iuj —>■ uj+i5) is the retarded self-energy 

matrix. We define the matrix spectral function, A(k, w), 
via 


G(k, iuj) 



dio 


,A(k^ 

iuj — uj' 


(13) 


such that 


A(k,a;) = ^(G«(k,a;)-G^(k,a;)) (14) 


where = G^^ is the advanced Green’s function. 
Since our calculation will only require A(k, w —>■ 0), we 
need only calculate 


G^(k,0) 


iFo-l-^ —vk —Afe 
—vk^ zFo -I- p, 0 
-Afc 0 fFo-^ 

0 —Afc vk'^ 


0 ■ 

—Afc 
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(15) 

where k^ = kx i iky. Here we have taken a simple form 
for the zero-frequency self-energy matrix, E^(a; —>■ 0) = 
—fFol, where Fq is a scalar constant, the impurity scat¬ 
tering rate. In general, the full 4x4 self-energy matrix 
can be calculated for a particular disorder model, but 
this simple model captures the essential physics and es¬ 
tablishes an energy scale for disorder. Performing the 
inversion in Eq. m yields the zero-frequency, matrix 
spectral function 


A(k,0) 


(AqUo- + AiCTi + ^2(72) Ir 


A 


den 


(16) 


where is the intra-block (spin) 2x2 identity matrix. 
It- is the inter-block (particle-hole) 2x2 identity matrix, 
and 

Aq = To (Fq-|-- f -f Afc) (17) 

Ai = 2rofJ,vkx 
A 2 = 2roMvky 

Aden = TT (Fq -f (vk — ^)^ -f Afc) (Fq -|- ( — vk — /i)^ -f Afc) 


B. Thermal Current Operator 


G(k, iuj)-^ = G°(k, iuj)-^ - E(iuj) (10) 

such that 

G(k,ioj) = [icul — S(ioj) — Hk] ^ (H) 


To derive an expression for the thermal current density 
operator in this system, we generalize the approach devel¬ 
oped for the s-wave superconductor case by Ambegaokar 
and Griffin^ and adapted for the d-wave superconductor 
case in Ref. |4^ We begin by expressing the Hamiltonian 
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in terms of the coordinate-space field operators, 
and V' 4 .(x), such that 

H = Ho + Hi (18) 

Ho = J (fx {-iva • V - /i) ^ 

-ffl = ^ Jd^X yd2y^t^V'y;3^(x-y)V'y/3V'xa 

where a and /3 are spin indices over which summation is 
implied, y(x —y) is the effective potential that gives rise 
to the proximity-induced superconductivity, and we have 
adopted a compact notation whereby 'i/'a = = V'c((x) 

and ij^yp = '4’pi'y)- Performing the matrix multiplica¬ 
tions, Hq takes the form 


Ipl = -'CCl+V't + Wxi’i 
■01 = 

0j = -vd~4i\. - iii\ipx (21) 

where we have defined 

(fx ='f{^) = J d^rV{r - (22) 

The thermal current density operator, j'‘(x), is ob¬ 
tained via continuity with the thermal density operator, 
/i(x). 


Ho 








iv (d -b 5 + 0101 ) “ MV'a0a (19) 


where 5+ = and the second equality is the result 

of integration by parts. Equations of motion for the field 
operators are obtained by noting that 


/i(x) =-V • j'‘(x) (23) 


Since we have written our Hamiltonian such that all ener¬ 
gies are measured with respect to the chemical potential, 
/i(x) is equal to the Hamiltonian density operator and 
therefore defined via 


i'ipa = [i’a,H] iipl = [ipl, H] (20) 

and applying fermion anticommutation relations. Doing 
so, we find that 



(24) 


= -vd 01 -b l(b>a; 0 t 


and expressed as 


J 


/i(x) = - 


(0j5 04,-5 0|04 -b 0|5+04- - 5+0J0-4) - M0l 0a + \Jd^y (25) 

where we have taken Ho to be the average of the first and second lines of Eq. m- Taking the time derivative 
breaking the result into two pieces, we write 

(i(x) = Fa + Fb 

where 

Fa = (0|5“04, - 5“0)f04. -b 0|5+0-i- - 5+0|0t -b 0|5“04 - 5“0|0^ -b 0|5+0-4 - 5+0|0- 


and 


(26) 


Fb = - (fyV{y-x) U 


(^lV'y/3V’y/30a + 0l0y/3V'y/30a + 0l0y/3^?//3^a + 


- M0i0a - M0i0a 

(27) 

(28) 

The first piece. Fa, can be reorganized by using the equations of motion (1^ to sub in for the dotted field operators, 
regrouping terms, and then applying the equations of motion again. Doing so, we find that 

+ 9+(0|0t) - 5"(0l04,) - 5+(0|0f) - J(fyV{y- x) (0l0j^0y/30a + 0l0j^0y/30a) (29) 

Combining this with Fb and applying the continuity equation (1231) . we see that it is natural to write the thermal 
current density operator as the sum of two terms 

j'" = ui -b U2 (30) 


Fa =— 

2 L 
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where 


V • Ui = — 


IV 


2 L 


d ('0IV'4,)+ 9+('!/’j'0t) - ^ (V’t'04.) - 


(31) 


V • U 2 = 


<Py ^(y - x) + V'lV’y/jV’y/sV’a) (32) 


Expansion of the operators reveals that the right-hand-side of Eq. dSD is easily expressed as a divergence. Doing 
so, we extract 


Ui = - 


IV 


2 1 


'I'Pi + X - i (i/'jV’t - i’l 


^ - i’l 


(33) 


which, in 4 X 4 Nambu notation, becomes 


ui(x,t) = 


4'^(7T34' — dt^CTTa'i' 


(34) 


where = 'P^(x, t) = ['0|, —V'tl ^ ~ + ^‘^V- Fourier transforming in space and time yields 

Ul(q,f2) = (‘^ + §) 


(35) 


where we have used the shorthand = 4'(k, uj) and '^k+q = (k -f q, w -b D). 

To obtain U 2 , we take the space-time Fourier transform of Eq. (ESI). Doing so yields 

iq-U2(q,D) = \J(fxcfydtV{y-x){e~^^'^-e~^^'y) + V’L^y/3^y/3^a;a) = Xi-bX2-Yi-l2 (36) 

where we have labeled each of the four resulting terms: Xi, X 2 , Yi, and ^ 2 - Inserting a Fourier representation for 
the potential and each of the field operators, the Xi term takes the form 

Xi = ^ ^ ^ i^iVk^cl^^cl^pCkspCkia S{k4 - ki-k^- q)5{k3 - k 2 + k5)S{uJi + UJ 2 - UJ 3 - UJ 4 + n) (37) 

All.. .fcs ■<■^4 


Making a mean field approximation, retaining only the 
terms for which the average values are over (k f, —k 4 .) 
pairs (reduced approximation), and noting that 
is an even function of w, this becomes 

X4=^J2ico-mlcL,■^cl,^ (38) 

kuj 

where 

Afc = ^ Ffc_fc'(c^,.|-cl^.,j^) (39) 

k'uj' 


In the g —>■ 0 limit, 

dXk 

Xk+q - A/c = q • = q • VAfe (41) 

where VAfc is the slope of the order parameter (the gap 
velocity) at k. Plugging into Eq. (HOI) and taking the 
D —>■ 0 limit, we find that 

U2(0,0) = ^ (w + ^ j VAfc + C-klCk+qt) 

kio ^ ^ 

(42) 

which, in the 4x4 Nambu notation, becomes 


is the superconducting order parameter for the interface 
state. Repeating this calculation for X 2 , Yi, and Y 2 , and 
taking A^ to be real, we find that 


q • U2(q, D) = ^(Afc+, - Afc) 


kcj 


+ (w + n)c_klCk+<rt\ (40) 


L12V0,0; - 2 

kuj 


Thus, in the g, D —>■ 0 limit (which is the limit where we 
will need it), the thermal current density operator is 


r(0,0) = i^4/l[a;-b^)vM4'fc+, (44) 

ku) ^ ^ 
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k, co + n 



k,(M 


trace is over Nambu space, and the factor of 1/2 out front 
compensates for the particle-hole double counting that is 
inherent in our 4x4 Nambu formalism. Inserting a ma¬ 
trix spectral representation, as defined in Eqs. (HI and 
da, for each of the Green’s functions, this becomes 


FIG. 1: Feynman diagram representing the bare bubble ther¬ 
mal current-current correlation function, tlre(*II). On each 
vertex sits a thermal current density operator, j"’. Each prop¬ 
agator line denotes a Green’s function dressed with disorder 
self-energy, G(k,iLj) and G{k,iuj + iQ). 


where 


vm = vaT3 + VAk Ti (45) 

is a vector in coordinate space and a matrix in our 4x4 
Nambu space. Here, the first term derives from the mass¬ 
less Dirac spectrum of the TI surface and has inherited 
the interesting spin structure thereof, while the second 
term derives from the the d-wave order parameter of the 
proximity-induced superconductivity. 


C. Thermal Conductivity 

With the spectral function and thermal current density 
operator in hand, we can proceed to calculate the ther¬ 
mal conductivity in the zero-temperature, zero-frequency 
limit. For d-wave superconductors, this limit is known as 
the universal limit because thermal conductivity has been 
shown to be insensitive to disorder in this regime 42r— 
We can calculate the thermal conductivity tensor, k(T), 
for the case at hand by appealing to the fluctuation- 
dissipation theorem as expressed in the Kubo formula^ 

M (46) 

T n^o T^n ^ ^ 

where we obtain the retarded current-current correlation 
function via analytic continuation from the Matsubara 
function. 

n^(D) (47) 

For simplicity, we proceed by calculating the bare bub¬ 
ble Feynman diagram shown in Fig.[Tl noting that vertex 
corrections have been shown to be small for the d-wave 
superconductor case^ and deferring to future work their 
calculation for the case at hand. Doing so, the Mat¬ 
subara thermal current-current correlation function takes 
the form 

iu) k 

X Tr [G(k, fa;)vMG(k, iw-I-*D) vm] (48) 


n,(zD) = 


X 


- f duji du!2 S{iV,) 

k 

Tr [44(k, a;i)vM^(k, a; 2 )vM] 


(49) 


where 


S{in) 


I X—^, zD. 9 I I 

-> (iwH- Y - 

B 2 zu; — wi zw -b zD — 

lUJ 


(50) 


Evaluating the Matsubara sum via contour integration 
(see Refs. and IU for a discussion of the technical 
points) and continuing zD —> D -|- zd, we obtain the re¬ 
tarded function 


SR{n) = 


{uji -b D/2)^nF(wi) — (a;2 — VL/2YnF{uJ2) 
UJi — UJ2 -b D -b zd 


(51) 

where npiuj) = -b 1) is the Fermi function. Since 

the retarded and advanced Green’s functions are hermi- 
tian conjugates, the spectral function defined in Eq. da 
must be hermitian 

t t r^R r^R r^A 

. 4 . Lt — Lx Lt — Lt Lt — Lt . 

A' = -z-= -i -= i -= A 


27T 


27T 


27T 


(52) 

And since vm is also hermitian, the trace in Eq. (PI) 
must be real 


Tr [AiVmA 2 Vm]* = Tr (AiVmA 2 Vm) 
= Tr 


(AiVmA2Vm)^ 

= Tr 


VmA2VmAi] = 

Tr [AiVmA 2 Vm] 


(53) 


Therefore 


Imn^(D) = f dwi di02 ImS'_R(D) 


X Tr [A(k,a;i)vMA(k,a; 2 )vM] (54) 


where 




ImS'fl(D) = 7r(a;i-b—) (rz_F(a;i-bD)-nF(wi)) 

X Siuji -b D — a;2) 


(55) 


Plugging into Eq. (HSl) and taking the D —^ 0 limit yields 
an expression for the thermal conductivity tensor. 



^Tr [A(k,a;)vMA(k, w)vm] (56) 

k 


where /3 = l/fc^T, the w-sum is over fermionic Matsub¬ 
ara frequencies, the fc-sum is over the Brillouin zone, the 


X 
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In the zero temperature limit, {uj/TY{—dnF/dui) is 
sharply peaked at w = 0. Thus, evaluating the integral 



we find that the universal limit thermal conductivity ten¬ 
sor takes the form 

^ Tr [A(k, 0)vM^(k, 0)vm] 

k 

(58) 

where A(k, 0) is the spectral function that we evaluated 
in Eqs. dill) and (13. 

Introducing the shorthand TrAvAv for the trace in the 
above expression and plugging in for A(k, 0) via Eq. (fTHll 
and for vm via Eq. (gSI) yields 


Ko _ k{T) 
T - T 


TyAvAv = 




Tr 


(Nir {vaTs + VAfc IcrTl))^ 


= ^^Tr -I-VAfeVAfciV^l^] (59) 

^den 


where N = Aola + Aiai+A 2 (J 2 and we have made use of 
the multiplicative properties of the particle-hole (r) Pauli 
matrices. Noting that a = aix + a 2 y, making use of the 
multiplicative properties of the spin (a) Pauli matrices, 
evaluating the trace, and plugging back into Eq. (1^ . we 
find that 


in for Aq, Ai, A 2 , and Aden from Eq. (HU, and restoring 
h in the prefactor, we obtain the following expression for 
the thermal conductivity in the zero temperature limit 


Ko k% 3 


VAfc VAfc (Pk + Qk) 

k k 


(63) 

where Pk = Al/A\^^ and Qk = (Af -f Al)/A\^^ take the 
form 


Pk = 


ho/Tr 


To/tt 


4 [rg-f (ufc - ^)2 + A2 Tl + {-vk- nf a 

(64) 


_ ho/Tr _ rp/TT _1 

Pp -I- {vk - nY -f Pq -I- {-vk - yY -I- A^ 

(65) 

Note that this result depends on integrals of the squares 
of sums and differences of Lorentzians centered about the 
zeros of the two branches of the quasiparticle excitation 
spectrum, Eq. (0), of width given by the impurity scat¬ 
tering rate. For /r ^ Pp, /Ip is dominated by impurity- 
induced quasiparticles in the vicinity of the zeros of the 
(-I-) branch. For /i <C —Tp, the (—) branch dominates. 
For 1^1 <C Tp, both branches contribute. 



IV. ANALYTICAL RESULTS 
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\xx + yy)'^ 


42 

k '^den 


■{xx-yy)^- 


42 

k ■^den 


\xy + yx)^ 


+ VAfcVAfc 


2^1 ^2 
“42 

k '^den 

AI + AI + AI 
Ain 


In both the large-|/r| and small-|/r| limits (|^| Pp and 
1^1 <C Pp) the quasiparticle excitation spectrum simpli¬ 
fies, as described in Sec. IIIBl and can be linearized about 
nodal points in fc-space. As a result, in these limits, we 
can obtain simple, closed-form expressions for the zero- 
temperature thermal conductivity. This is shown in the 
following sections. 


(60) 


A. Large-1/rI Limit 


Since Ak is of (1^2 _y 2 symmetry, it must be an even func¬ 
tion of both kx and ky. Therefore, as defined in Eq. (HU, 
Ap and Apen are even functions of kx and ky while Ai is 
odd in kx but even in ky and A 2 is even in kx but odd in 
ky. As a result 


\ "v 2 A 1 A 2 f dkx f dky 2 A 1 A 2 

k ^den J 2tT J 2tT A^^^ 


And since exchange of kx for ky sends A^, to —Ak, it 
leaves Apen invariant but exchanges Ai for A 2 . Therefore 


Al 

2-^ a‘^ 2 ^ ^2 

k '^den ^ "^dei 


(62) 


Thus, only the first and fourth terms in Eq. (16011 survive. 
Noting that xx+yy is just the identity tensor, 1, plugging 


For /r ^ Tp, the Lorentzians in Eqs. (|M|) and (113 
are sharply peaked about four nodal points, located at 
±kx = Aky = /i/-\/2u, and are well separated from each 
other. We can therefore replace the fc-sum in Eq. ((63ll by 
the sum of four integrals over local scaled coordinates, pi 
and p 2 , defined about the nodal points 


4 

E-E 


fc i=i 


d^k A 

^ h 


d?p 

{2'kYvv/^ 


( 66 ) 


where the integrals can be extended to infinity because 
the integrands are so sharply peaked about each node. 
Here, pi = vki and p 2 = VAk 2 , and at each node, fci 
and ^2 point, respectively, perpendicular to and parallel 
to the local Fermi surface, with ^2 in the direction of 
increasing A^. In terms of these scaled coordinates, A^ r:! 
P 2 and vk = p + pi, so vk — p = Pi and —vk — p = 




































— {2fi + pi) ~ —2p. Therefore, since p ^ Tq, the second 
Lorentzian can be neglected with respect to the first in 
both Eq. (l64ll and Eq. (l65ll and we find that 


Pk 



( To/tt 


2 


(67) 


where p = yjp\ + p\- Evaluating the integral 


/ 


(Pp 


( rp/TT 

\rl + p^ 


2 


1 

47r3 


( 68 ) 


and noting that the sum over nodes of the outer product 
of VAfe with itself at each node is 


E' 


= 2vll 


(69) 


we find that 




1 

VVA 


1 1 
4 47r3 


1 

4tt^vva 


(70) 


E Q O 1 11 "U A 1 

VAfcVAfc {Pk + Qk) = 27;|1- 2’-’—^ = —- 

VVA 4 47r'^ 47r'^z;x;A 

k 

(71) 

Therefore, the thermal conductivity tensor reduces to a 
scalar, )^q = kqI, with the simple form 


T ~ 2 3n \VA ^ V J 


(72) 


The same result is obtained for p <C —Tp, where it is 
the Hrst Lorentzian in Eqs. ([M]) and (Ell) that can be 
neglected with respect to the second. Note that this ex¬ 
pression is independent of disorder and is only a function 
of the velocity anisotropy, v/va, which depends on both 
p and material parameters. Note also that this is exactly 
half the value obtained (per layer) for the case of an or¬ 
dinary d-wave superconductor 4^ This is because, unlike 
the d-wave superconductor case where the electron dis¬ 
persion is spin-degenerate, here the TI surface state is 
nondegenerate and only one of the two branches of the 
quasiparticle excitation spectrum (Eq. (El) contributes 
to the thermal conductivity. For ^ Tp, the (-I-) branch 
(first Lorentzian) contributes. For p <C —Tp, the (—) 
branch (second Lorentzian) contributes. This factor of 
two is a clear and measurable demonstration of the sense 
in which the TI surface topological metal is “half” of an 
ordinary 2D electron gasi^. 


B. Small-|/i| Limit 

For \p\ <C Tp, the four anisotropic nodes of the prior 
section have coalesced into a single isotropic node at the 
origin of fc-space. The first and second Lorentzians in 


Eqs. El and (1651) are approximately equal and peaked 
at the origin. The fc-sum in Eq. (1631) can be replaced by 
a single integral about scaled coordinates, pi = vk^ and 
P2 = vky, and extended to infinity. 



In these scaled coordinates, vk = p = \Jp\ -\- p\, and 
since \p\ <C Tp, vk — p ^ p and —vk — pm —p. As 
long as Ajt vanishes fast enough with decreasing fc, as 
per condition (2) of Sec. IIIBl and 'VAk'^Ak can be 
neglected in Eqs. (I631I65I) compared to larger terms. As 
a result, the two Lorentzians add in Pk and cancel out in 

Qk- 


Once again making use of the integral in Eq. (1551) . we 
find that 


V p = — /" ( JEaV = 

^ v'^ J (27r)2 


(75) 


VAfcVAfc (Pk + Qk) ~ 0 (76) 

k 


Therefore, the thermal conductivity tensor again reduces 
to a scalar, now with an even simpler form 


12 - M.i 

T ~ 3h 2 


(77) 


Here, both branches of the quasiparticle spectrum have 
contributed to the thermal conductivity, and one obtains 
precisely the result one would expect for a single isotropic 
massless Dirac node. This expression is clearly indepen¬ 
dent of disorder and is just the standard d-wave super¬ 
conductor result^ for an anisotropy ratio of one, divided 
by a factor of four since there is only one node here rather 
than four. 


V. NUMERICAL RESULTS 

We would now like to look beyond the large-1^| and 
small-1/r I limits and consider the transition between them 
by numerically evaluating Eqs. (I63II65|) as a function of 
p. This is easily done, but unlike the large and small |/r| 
limit calculations which were model independent (aside 
from the two conditions in Sec. ttm, this calculation re¬ 
quires a model for A^, the proximity-induced supercon¬ 
ducting order parameter of the Tl-dSC interface state, 
and its results will necessarily depend (in the details) on 
that choice of model. Since we are primarily interested 
in understanding the essential physics of this transition, 
without delving too deeply into the material-dependent 
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details, we proceed by considering the following simple 
and rather standard expression for the order parameter 
of a generic d-wave superconductor 


A/e = (cos k^a — cos kyo) 


(78) 


which yields a gap velocity at k of the form 

^ _ aA/e Aoa - , • , 

VAfe = ^ (- sm kxa x + sm fcy a y) 


(79) 


Here, we have introduced two new model parameters, the 
gap maximum Ag and the lattice constant a. Expressing 
all lengths in units of a and all energies in units of n/a, we 
define dimensionless parameters fi = fj.a/v, Eg = Toa/v, 
and Ag = Aga/n, as well as a dimensionless wavevector 
with components Zi = kxa and Z 2 = kya. Doing so, plug¬ 
ging Eqs. (|781I79I) into Eqs. (I631I65I) . and noting that all 
terms not proportional to the identity tensor integrate to 
zero, we find that the universal-limit thermal conductiv¬ 
ity tensor reduces to a scalar and takes the convenient 
form 


Y 


k'^ 


2 , r 


d^z 


3fi J Stt 
A2 


{L{z)+Li-z)Y 


sin^ Zi {L{zY + E(—z)^) 


(80) 


where 


and 


L{z) = 


Tl + {z-Yy+A{zY 


A(z) = (coszi — COSZ 2 ) 


(81) 


(82) 


The fc-space integral is easily computed to obtain kq/T 
as a function of fi. Results for Ag = O.lv/a and Tg = 
Q.Q\v/a are plotted in Fig. [5] alongside the large and small 
1^1 limits. (For the large-1^| plot, we have used the model 
introduced in Eqs. (I781I79I) to obtain the nodal anisotropy 
ratio as a function of /i, v/vi\ = [(Ag/-\/2) sin(/i/-\/2)]“^, 
and used that as input to Eq. dll-) Our numerical re¬ 
sult matches the large-^ expression for |/r| Tg, peak¬ 
ing with decreasing |^|, before plunging down toward the 
small-|^| value for |/r| <C Eg. 

This behavior is best understood by considering the 
evolution of the fc-space structure of the integrand of 
Eq. (1501) as a function of //, as shown for a series of fi 
values in Figs. |3] and 01 The upper panel of Fig. |3] il¬ 
lustrates the structure of the large-/x limit. Here, for 
fi = 7r/\/2, the integrand is peaked within Eg of four, 
well-separated, anisotropic nodal points. Equal-intensity 
contours are (nearly) elliptical, squeezed in the direction 
parallel to the local Fermi surface. In the middle panel, 
fi is reduced by a factor of 2, which draws the nodes 
closer to the origin. The peaks are still well-separated. 



FIG. 2: Calculated universal-limit thermal conductivity, 
Ko/r, as a function of chemical potential, fi. Solid curve de¬ 
notes numerical solution of Eqs. (I80I82I) for parameter values 
Ao = O.lu/o and Eq = O.Olu/a. Solution matches our large- 
IpI expression (dashed) for |p| ^ Eq, then reaches a maximum 
at an intermediate value of |/r| before decreasing toward the 
value of our small-|/i| expression (dotted) for |p| <C Eq. 


but less so than before, since the radius of the Fermi 
circle has decreased and the nodal anisotropy ratio has 
increased. Thus, the peaks have begun to curve around 
the Fermi circle, toward each other, and the independent 
node approximation used to derive the large-1^| expres¬ 
sion of Eq. (1721) has begun to break down. In the lower 
panel, fi is reduced by an additional factor of 10. Now the 
independent node approximation has completely broken 
down, and the four anisotropic peaks have curved into 
each other, forming an annulus of width Eg about the 
Fermi circle. With decreasing /i, the radius of this annu¬ 
lar peak decreases, resulting in the decrease of Kg/T seen 
in Fig. 01 We reproduce this image, zoomed-in about the 
annulus, in the upper panel of Fig. [3] In the middle panel 
of that figure, fi is reduced by another factor of 10, such 
that it is nearly equal to Eg. Now the width and radius 
of the annular peak are nearly equal. As fi decreases 
further, the system is tuned toward the Dirac point in¬ 
herited from the TI surface state and the isotropic node 
at the origin is revealed. For |/i| <C Eg, the annular peak 
blurs into a single isotropic peak at the origin, of width 
Eg. This is shown in the lower panel where fi = 0. The 
integral over this single isotropic peak recovers the small- 
\fi\ value of Eq. (EZl). As fi becomes negative, the process 
reverses, dominated now by the (—) branch of the quasi¬ 
particle excitation spectrum instead of the (-)-) branch. 
All else is the same, so kq/T is even in fi. 

Results for five different values of the impurity scat¬ 
tering rate. Eg, are shown in Fig. [SI Note that in 
both the large-1^1 and small-|/r| limits, kq/T is disorder- 
independent. The transition between these limits does, 
however, depend on disorder, with the peaks of the kq/T 
vs fi curves smoothed out for greater disorder. This ef- 
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FIG. 3: Evolution of the fc-space structure of the Ko/r inte¬ 
grand. (a) [fi = Large-I^l limit. Four, well-separated, 

elliptical peaks within Fo of the nodal points, (b) [/i = 

Nodal peaks closer to the origin, more anisotropic, and curv¬ 
ing around the Fermi circle, (c) [fi = Nodal peaks 

have merged into an annular peak of width Fq. 


FIG. 4: Further evolution of the fc-space structure of the ko /T 
integrand, zoomed-in by a factor of 20. (a) [/r = 

Closeup view of the same annular peak shown in Fig. I^c). 
(b) [/r = Width and radius of the annular peak 

now nearly equal, (c) [/r = 0] Small-1^| limit. Annular peak 
blurred into single, isotropic peak within Fq of the origin. 
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FIG. 5: Disorder dependence of calculated universal-limit 
thermal conductivity, ko/T, as a function of chemical po¬ 
tential, fi. We plot numerical solutions of Eqs. (I80II82I) for 
Ao = Q.lv/a and five values of the impurity scattering rate, 
Fq. Results are disorder-independent in both the large-|p| 
and small-1/i I limits. The transition between limits depends 
on disorder, with the peaks more prominent for smaller Fo, 
smoothing out with increasing disorder. 


feet can be understood in terms of our integrand analysis 
(above). As |/r| decreases from its largest values, increas¬ 
ing anisotropy ratio yields increasing kq/T via our large¬ 
ly] expression, Eq. (17^ . But for greater disorder, the in¬ 
dependent node approximation that defines the large-1^| 
limit breaks down sooner, as the four anisotropic peaks 
broaden with growing disorder and merge together ear¬ 
lier, limiting the enhancement of nojT with increasing 
anisotropy ratio. The resulting annular peak is of greater 
width for greater disorder and therefore blurs into a single 
peak sooner, ushering in the small- |/i| limit as its radius 
becomes smaller than its width. 


VI. CONCLUSIONS 


In this paper, we have calculated the universal-limit 
thermal conductivity, kq, as a function of chemical poten¬ 
tial, fjL, due to quasiparticle excitations of the proximity- 
induced superconducting state at the 2D interface of a 
topological insulator and a d-wave superconductor. In 
both the large-|/r| and small-|/r| limits, we have obtained 
simple closed-form expressions for kq/T, combined here 
from Eqs. (1771) and ((7^ 


= + l/^l » To 

^ I 1 for |/r| <C Eq 


(83) 


where v is the slope of the isotropic Dirac cone inher¬ 
ited from the TI surface state, v/va is the ^-dependent 
anisotropy ratio of the four anisotropic Dirac cones of the 
proximity-induced d-wave superconducting state, and Eo 


is the impurity scattering rate, the energy scale char¬ 
acterizing disorder in the system. Note that the large¬ 
ly! expression is exactly half the value obtained^ (per 
layer) for an ordinary d-wave superconductor: = 

{kg/3h){vF/vA+VA/vF)- This is an overt demonstration 
of the sense in which the underlying topological metal 
is “half” of an ordinary metali^, and comes about be¬ 
cause, for large |^|, only one of the two branches (positive 
or negative) of the isotropic Dirac cone contributes at a 
time. Eor |/i| <C Eq, both branches contribute, but the 
four nodes have coalesced into one isotropic node at the 
origin of A:-space. Thus, the small-|p,| expression is equal 
to the standard dSC value (with anisotropy ratio equal 
to one), divided by four (since there is only one node in¬ 
stead of the usual four): (1 -|- l)/4 = 1/2. While kq/T is 
disorder-independent in both of these limits, the transi¬ 
tion between them, as a function of /r, depends on disor¬ 
der. And furthermore, it depends, in the details, on the 
functional form of the proximity-induced order parame¬ 
ter, Afe. Adopting a simple model for (Eq. (1751) '). we 
have calculated kq/T across the full range of /r, for differ¬ 
ent levels of disorder, as shown in Eig. [S] As /i decreases 
from its maximum value, the four nodal peaks of the in¬ 
tegrand in Eq. (1801) become more anisotropic, resulting 
in an increase in kq/T, as per our large-|/i| expression. 
But they also move closer together, eventually merging 
into an annular peak about the Eermi circle. Along the 
way, the independent node approximation that defined 
the large-|/r| limit breaks down, and kq/T reaches its 
maximum value, decreasing as /r decreases further and 
the Eermi circle shrinks. Einally, as fj, gets smaller than 
Eq, the annular peak blurs into an isotropic nodal peak 
at the origin, and kq/T reaches its minimum at the value 
given by our small-|^| expression. As shown in Fig. [51 
the peaks in the kq/T vs ^ curve are more pronounced 
for smaller Eg, smoothing out with increasing disorder. 

Note that we have assumed herein that the bulk band 
gap of the topological insulator extends well above and 
below the Dirac point of the surface state, such that 
could be varied over a wide range of energies without 
accessing the bulk valence or conduction bands. In real 
materials, the available energy windows may be more re¬ 
stricted. We have also assumed that the chemical po¬ 
tential can be accurately controlled, via gating, doping, 
or other means, and that proper contact can be made 
to the Tl-dSC interface. Both may present experimental 
challenges. 

Our focus in this work has been on the evolution with 
changing chemical potential of the massless Dirac quasi¬ 
particle excitations of the Tl-dSC interface state. Results 
shed light on the essential features of low-temperature 
thermal transport due to these quasiparticles. Further 
theoretical development, including incorporation of a 
subdominant spin-triplet order parameter, a more real¬ 
istic disorder model, and vertex corrections to our dia¬ 
grammatic calculation, are left for future work. 
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